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I. INTRODUCTION 



A. PURPOSE 

This thesis describes the development of a microcomputer 
oriented program called MSCOP (Microcomputer Software for 
Constrained Optimization Problems) for constrained optimiza- 
tion of engineering design problems. Problems which can be 
solved by the MSCOP are nonlinear programming problems 
arising in several areas of machine and structural design, 
such as the minimum weight design of structures subject to 
stress and displacement constraints [Ref. 1]. 

In recent years, several powerful general purpose opti- 
mization programs have become available for engineering 
design problems, e.g., COPES/CONMIN [Ref. 2], and ADS-1 
[Ref. 3], These programs can handle a wide range of design 
problems and contain a variety of solution technigues. 
Also, several programs are available that include optimiza- 
tion in an integrated analysis / design code, e.g., ACCESS, 
ASOP, EAI, PARS, SAVES, SPAR, STARS and TSO [Ref. 4]. All 
of the above optimization programs are written in FORTRAN, 
and are built for use on a mainframe computer. Their use can 
be cumbersome, especially for the occasional user. Since 
many engineers are now using microcomputers, there is a need 
to develop an optimization program contained in a microcom- 
puter software package for use on microcomputers. This 
thesis fills that need by developing a compact program 
written in a standard EASIC language suitable for a wide 
range of microcomputers. 
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B. I HP LEHS STATION 



The nature of an optimization program depends on the 
computer and programming method available. The MSCOP soft- 
ware is designed for use on a microcomputer. However, for 
the speed of development and testing, MSCOP was developed on 
the IEM 3033 computer at the . R. Church Computer Center in 
Naval Postgraduate School, and was written in WEASIC 
(Waterloo Basic) Version 2.0. 

To make sure that the program is easily portable to a 
microcomputer, only standard BASIC commands and functions 
are used. For example, FOR I = 1 TO MDB ... NEXT I, G0SU3 
etc., were used. The commands and functions not available 
in all variations of EASIC are avoided, for example, TRN(A), 
MAT (A), etc. 

MSCOP provides design engineers with a convenient tool 
for optimization of engineering design problems with up to 
20 bounded design variables and as many as 50 ineguality 
constraints. 



C. GENERAL OPTI MI Z ATICN MODEL 



The general optimization problem to be solved is of the 
form : Find the set of design variables _X that will 



Mini mi ze 
Subject to 



F(X) 

G (X) < 0 

j 

1 u 

X < X < X 
i - i - i 



^ _ 1,..., m 



i = , n 



( 1 . 1 ) 
( 1 . 2 ) 

(1.3) 



where X is referred to as the vector of design 
F (X) is the objective function which is to be 
G (X) are ineguality constraint functions, and 
are lower and upper bounds, respectively, on 



variables. 

minimized . 

I , a 

X t and X c 

the design 



variables. Although these bounds or "side 
could be included in the inequality constraint 
Eg (1.2) , it is convenient to treat them separa 
of their special structure. The objective f 
constraint functions nav be nonlinear, explicit 
in X. However, they must be continuous and 
continuous first derivatives. 

In general engineering optimization problems 
tive to he minimized is usually the weight or 
structure being designed while the constraints 
cn compressive stress, tensile stress, Eule 
displacement, frequencies (eigenvalues) , etc. 
p.264 ]. Equality constraints are not included b 
inclusion complicates the solution techniques an 
engineering situations, equality constraints are 

Most optimization algorithms require that 
value of design variables X° be specified. Be 
these starting values, the design is iterative 
The iterative procedure is given by 



q+1 q q 

X = X + a* S 



constraints" 
set given by 
telv because 
unction and 
or implicit 
should have 

, the objec- 
volume of a 
gives limits 
r buckling, 
[Ref. 5 : 
ecause their 
d because in 
rare. 

an initial 
ginning from 
ly improved. 

(1.4) 



where q is the iteration number, S is a search direction 
vector in the design space, and a* is a scalar parameter 
which defines the amount of change in X. At iteration q, it 
is desirable to determine a direction S which will reduce 
the objective function (usable direction) without violating 
the constraints (feasible direction). After determining the 
search direction, the design variables, X, are updated by Eq 
(1.4) so that the minimum objective value is found in this 
direction. [Ref. 6]. 

Thus, it is seen that nonlinear optimization algorithms 
for the general optimization problem based on Eg (1.4) can be 
separated into two parts, determination of search direction 
and determination of scalar parameter a*. 



1 2 



D. 0 EG ANIZATION OF IBIS THESIS 



This chapter has stated the purpose of the thesis and 
has put the general concept of engineering optimization into 
a preliminary perspective. Chapter 2 will describe the 
essential aspects of the optimization algorithm used in 
MSCOP such as finding a search direction, the one- 
dimensional search and convergence criteria. Chapter 3 
describes program usage. In chapter 4, there are three 
examples which are solved by the MSCOP. Summary and conclu- 
sions are given in chapter 5. The program is listed in the 
appen dix . 
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II. OPTIMIZATION ALGORITHM 



A. INTRODUCTION 

There are many optimization algorithms for constrained 
nonlinear problems such as generalized reduced gradient 
method, feasible direction method, penalty function methods. 
Augmented Lagrangian multiplier method, and sequential 
linear programming. The feasible direction method is chosen 
for development in this thesis for three main reasons. 
First it progresses rapidly to a near optimum design. 
Second it only requires gradients of objective and 
constraint functions that are active at any given point in 
the optimization process [Bef. 7]. Third, because it main- 
tains a feasible design, engineer cannot fail to meet safety 
requirements as defined by the contraints. However, the 
method does have several disadvantages in that it is prone 
to "zig-zag" between constraint boundaries and that it is 
usually does not achieve a precise optimum. This method 
solves the nonlinear programming problem by moving from a 
feasible point (can be initially infeasible) to another 
feasible point with an improved value of the objective 
value. 

The following strategy is typical of feasible direction 
method : Assuming that an initial feasible point X° is 

known, first find a usable-feasible direction S. The algo- 
rithm for this is similar to linear programming and comple- 
mentary pivoting algorithms. Having found the search 
direction, a move is made in this direction to update the X 
vector according to Eg (1.4) . The scalar a* is found by a 
one-diraensicnal search to reduce the objective function as 
much as possible subject to constraints. That is MIN 
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( STARTS 







Figure 2.1 Algorithm for the Feasible Direction flethod. 



F(X+a*S) subject to G(X+a*S) £ 0. It is assumed that the 
initial design X° is feasible, but if it is not, a search 
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direction is found which will direct the de 
feasible region. After updating the X° vector, 
gence test must be performed in the iterative al 
convergence criteria used in this is impleme 
described in section C. The general algorithm u 
is given in Figure 2.1 

B. SEARCH DIRECTION 

In the feasible direction algorithm, a usabl 
search direction S is found which will reduc 
tive function without violating any constrain 
finite move. It is assumed that at any point i 
space (at any X) the value of the objective an 
functions as well as the gradients of these fu 
respect to the design variables can be calcula 
these gradients cannot usually be calculated a 
the finite difference method Eg (2.1) is used in 

ar (X) f (x+ ee ) - f (X) 

i 

ax e 

i 

where e is the ith unit vector 
i 

£ is a small scalar. 

In MSCOP , £ is 0.1% of the ith design v ar 

In the feasible direction algorithm, there 
one or more "active" constraints. A constraint 
"active" at X if g (X) Si 0. As shown in Figure 
constraints are active the standard steepest de 
tion S = - VF is used. 



sign to the 
the conver- 
gorithm. A 
ntation are 
sed in "SCO? 



e - feasible 
e the objec- 
ts for some 
n the design 
d constraint 
notions with 
ted. Since 
nalvtically , 
MS CO F . 



( 2 . 1 ) 



iable 






are 


usua 


11 y 


S(X) 


< 0 


is 


2. 1 , 


if 


no 


scent 


dir 


ec- 
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1 . Usa ble- Feasible Dire ctio n 



*2 




Figure 2.2 Osable-Feasible Direction. 

Assume there are NAC active constraints at X.. The direction 
S. is "usable" if it reduces the objective function, i.e., 

VF • S < 0 (2.2) 

Similarly the direction is feasible if for a small movement 
in this direction, no constraint will be violated, i.e., 

7G-S < 0 (2.3) 

This is shown geometrically in Figure 2.2 
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2 . 



Active Constrain ts 

It is necessary to determine if a constraint is 
active or violated in the feasible direction algorithm. A 
constraint G (X) <0 is "active" at X° if G (X°) *0. In order 

to avoid the zigzagging effect between one or more 
constraint boundaries, a tolerance band about zero is used 
for determining whether or not a constraint is active. From 

the engineering point of view, a constraint G (X) <0 is 

active near the boundary G (X) = 0 whenever ACC < G (X) < 7CC. 

ACC is the active constraint criterion and VCC is the 

violated constraint criterion in MSCOP . Assuming the 

feasible constraints are normalized so that G (X) ranges 
between -1 and 0 for reasonable values of X, the constraint 
G(X) <0 is considered active if G (X) > -0.1. The 

constraint is considered to be violated if G(X) > 0 .004. 
This is an algorithmic trick which improves efficiency and 
reliability of the algorithm. However, since in the one - 
dimensional search, all interpolations for constraint G (X) 
are done for zeros of a linear or quadratic approximation to 
G (X) in order to find a*, at the optimum the value of active 
constraints are very near zero, but may be as large as 0.004 
[Ref. 6]. From an engineering point of view, a 0.4 % 

constraint violation is considered to be acceptable. 



3 . Subopt i mizat ion Pro blem and Push^Of f Factors 

Zoutendijk [Ref. 8] has shown that a usable 
feasible direction S may be found as follows : 

Maximize ^ (2. 4) 

Subject to ; 

gF (X) - S + £ < 0 (2. 5) 
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j J 



( 2 . 6 ) 



7G(X)-S + 6 ^ < 0 

S bourded 



(2.7) 



Where scalar § is a measure of the satisfaction of the 
usability and feasibility requirements. The scalar 0j in 
Eg (2. 6) is referred to as the "push-off" factor which effec- 
tively pushes the search direction away from the active 




Figure 2.3 Push-Off Factor and Bounding of the S-Vector. 

constraints. In Eg (2. 6) , if the push-off factor is zero, 
the search direction is tangent to the active constraints, 
and if it is infinite, then the search direction is tangent 
to the objective function. It has been found that a 
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push-off factor is defined as follows gives good results 
[Ref . 5 : p. 167 ] : 



8 = 

j 



where 0. 



G . (X) 
3 

- ACC 




1 . 



( 2 . 3 ) 



To avoid an unbounded solution when seeking a usable 
- feasible direction it is necessary to impose bounds on the 
search direction S. Cne method of imposing bounds on search 
direction is to impose bounds on the components of S-vector 
cf form : 

- 1 < S. < 1 (2.9) 



This choice of bounding the S-vector actually biases the 
search direction. This is undesirable since we wish to use 
the push-off factors as our means of controlling the 
search direction. A method which avoids this bias in search 
direction is the circle as shown Figure 2.3 . The norm here 
is 

S • S < 1 (2.9.1) 

4 . S impl e S im plex -like Method for Search Direc tion 

Vanderplaats [Ref. 5: pp. 168-169] provides the 
matrix formulation which solves the above sub-optimization 
problem by using the Zoutendijk method. 



Maximize P-v 
Subject to ; 

A -y < 0 



( 2 . 10 ) 



( 2 . 11 ) 
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y-y < 1 

Where 



f s 

s 




0 


“1 






s 




0 


"2 






• 


P - 


» 


• 

s 




0 


~ n 






J 


• 


1 



( 2 . 12 ) 



(2. 13) 



V T G (X), 0 ' 

h 2 m. e 2 



A 



VG . (X), 
T 3 

VF (X), 




(2.14) 



and where j is the number of active constraints (NAC) 



When the solution to Eg (2. 10) through (2.12) is found, S may 
be normalized to some value other than unity, but the form 
of the normalization is the same. A solution to the above 
problem may be obtained by solving the following system 
derived from the Kuhn-Tucker conditions for that problem : 



[ 



B 




c 



u > 0 
i - 



v >0 u • v = 0 

i - ~ 



Wher e 



(2. 15) 



(2. 16) 
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(2. 17) 



B = -A* A 

I = Identity matrix (2.18) 

c = -A -P (2. 19) 

Atove system can be solved using a complimentary pivot algo- 
rithm. Choose an initial basic solution to Eq{2.15) is to 
be 



v = c, u = 0 ( 2. 20) 

where v is the set of basic variables and u is the set of 
nonbasic variables. If all v* t > 0, Eg (2. 16) is also satis- 
fied and problem is solved. If some v- L < 0, the solution 

procedure is as follows : 

Let be the diagonal element of the i-th nonbasic vari- 

able. 

1. Given the condition that some c is less then zero, 

we find max ( c ^ /3 ) which is the incoming row to the 

basis. 

2. The incoming column is changed to a basic 
column, the tableau is updated by a standard simplex 
pivot on B ^ • 

3. Until all c^> 0, repeat steps 1. and 2. 

4. When all c* >0, the iteration is complete. The 
value of u is now the desired solution. 

5. By using y = P~A T .u, we get the usa ble- feasible 
search direction S which is first NDV components of 

y- 

5 . In i ti a lly In feasible Desi gn s 

The method of feasible directions assumes that we 
begin with a feasible design and feasibility is maintained 
throughout the optimization process. If the initial design 
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is infeasible, then a search direction pointing toward the 
feasible region can be found by a simple modification to 
direction finding problem. 

A design situation can exist in which the violated 
constraints are strongly dependent on part of the design 
variables, while the objective function is primarily depen- 
dent on the other design variables. This suggests a method 
for finding a search direction which will* simultaneously 
minimize the objective while overcoming the constraint 
violations. These considerations lead to the following 
statement of the direction finding problem [Ref. 5 : 
pp. 171-172 ] : 



Maximize 



- vp (X) • S + ££ 



( 2 . 21 ) 



Subject to ; 



VG (X) • S + 0 £ <0 

j 



je j 



( 2 . 221 



S-S < 1 



(2.23) 



where J is the set of active and violated constraints, and 
where the scalar £ in Eg (2. 21) is a weighting factor deter- 
mining the relative importance of the objective and the 
constraints. Usually a value cf 5. > 10000 will ensure that 
the resulting S- vector will point toward the feasible 
region. Incorporating Eg (2. 21) and Eg(2.22) into the direc- 
tion finding algorithm reguires only that we modify the 
p-vector given in Eg (2. 24) and the A-matrix of Eg (2.25). 



P 



- VP(X) ' 

5 
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(2. 24) 



A 



2G (X), 0 

^ G 2 (X)/ 0 2 



V G.(X), 0. 

3 3 

0 <50 

3 “ 



(2.25) 



(2.26) 



We use the simple simplex-like method to find the 
search direction toward the feasible region. 



C. ONE-DIBENSIONAL SEARCH 

1 • M2 Violated Ccnstrai nts 

If no constraints are violated, we find the largest 
a* in Eg (1.4) from all possible values that will minimize 
the objective on S without violating any constraints, active 
or inactive. 

The procedure in MSCOE is as follows : 

1. Let aO, al, a2, a3 be the scalar in Eg (1.4) ccrre- 
sponding to pcints X0 , XJ, X2, X3, X4 . 

2. aC = 0 at given point X 0 . 

3. In order to get al, we can calculate the al to 
reduce the objective by at most 10% or to change each 
of the design variable X. by at most 10%. 

4. Update the design variables to XJ using Eg(1.4) . 

5. Evaluate the objective for XJ, and check the feasi- 
bility. If one or more constraints is violated, then 
al is reduced to al/2, and we go to step 4. 

6. In order to estimate a2, we can use the quadratic 
approximation with 2 points X, XJ. and the VF. 
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7. Update the design variables to X2 by Eg (1.4) and 
check the side constraints. 

8. Evaluate the objective and constraints. 

9. Now having 3 a's, and values of objectives and 

constraints for design variables XO, XJ , 12 are 

known, so by using 3-point quadratic approximation, a 
value of a3 is found. 

10. Update the new optimal point in search direction by 
Eg (1.4) . 

11. Evaluate the objective and constraints. 

12. Now choose last 3 values, al , a2, a3 and find a new 

a3 using 3-points Quadratic approximation 

13. Choose the a* among the 5 points which corresponds to 
the minimum objective function value with no-viclated 
ccnstrai n ts. 

2 • One or Mor e Con stra i nts V iol a te d 

If one or more constraints are initially violated, a 
modified usable-feasible direction is found. It is then 
necessary to find the scalar a* in Eg (1.4) which will mini- 
mize the maximum constraint violation, using the most 
violated constraint j, a good initial estimate for a* is 

-G . (X) 

a* = (2.27) 

VG (X) • S 

j 

Since the gradients of the violated constraints are 
known, the scalar which is required to obtain a feasible 
design with respect to violated constraint in the search 
direction, is given to a first approximation by Eq(2.27). 

The more detail procedure in ilSCOP is as follow ; 

1. Choose the most violated constraint j. 

2. Calculate a* for violated constraint j using 
Eg (2.27) . 
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3 . 



Update the design variables for a* and check the side 
constraints. 

4. If one or more violated constraints still exist, then 
calculate the derivative of objective, violated and 
active constraints and find a new search direction 
and then go tc step 1. Otherwise proceed with the 
optimization in the normal fashion. 

C. CCNVEBGENCE CRITERIA 

A desired property of an algorithm for solving a nonli- 
near problem is that it should generate a sequence of points 
converging to a global optimal point. In many cases, 
however, we may have to be satisfied with less faverable 
outcomes. In fact, as a result of non-convexity, problem 
size, and other difficulties, we may stop the iterative 
procedure if a point belongs to a described set, which is 
defined in MSCOP as fellows ; 

1. Q = [ X J |X° - X| < £ x |X°| } 

2. Q = { X 1 ! F (X°) - F (X) | < 4-|F(X0)| } 

2 ^ 

In MSCOP, the algorithm is terminated if a point X. is 
reached such that X & Q, f) Q x . £ x is 0.001 and is 

approximatly 0.001. Since in engineering design problems it 
is not necessary to find solutions with more than three 
significant digits. 
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III. HSCOP USAGE 



A. INTRODUCTION 

Since this MSCOP is written in WATERLOO BASIC Version 
2.0, it is verv convenient to use. The user must first 
formulate the design problem with the classical machine 
design criteria. Given the formulation of the design 
problem as a nonlinear program, the user then enters the 
problem as a part of a BASIC program. The user defines the 
objective function and constraint functions using EASIC 
statements. Other parameters are input as data : the number 
of design variables NDV, the number of inequality 
constraints NIQC, variable bounds an initial design value 
and a print control number. 

E. PROBLEM FORMULATION 

Generally, the experienced design engineer will be able 
to choose the appropriate objective for optimization 
depending on the requirements of the particular application. 
The physical phenomena of significance should first be 
summarized for the device to be designed. The appropriate 
objective can then be selected and constraints can be 
imposed cn the remaining phenomena to assure an acceptable 
design from all standpoints. However, the initial formula- 
tion for the optimization problem should not be more compli- 
cated then necessary and this often requires the making of 
some simplifying assumptions. [Ref. 9]. 

After completing the formulation of the design problem, 
the design engineer should be able to answer the following 
questions : 

1. What are the design variables ? 
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2 . 
3 . 
4. 



What is the objective function ? 

What are the inequality constraints ? 
What are the bounds on the variables ? 



The engineer is then ready to input the program to the 
MSCOP. However, additional study and preparation of the 
problem may be useful. In particular, redundant constraints 
should be avoided if possible. MSCOP will operate with 
redundant constraints but it will operate faster without 
them. Selection of an initial design point from which to 
start this program is important, since it affects perform- 
ance and running time. The user should use any available 
information which gives a good initial approximation. If 
side constraints exist, the user must be sure the initial 
values cf the design variables do not violate the side 
constraints. This program will automatically handle an 
initial design point which is infeasible with respect to the 
G (X) < 0 constraints. However, if the initial point does not 
violate these constraints, convergence will likely be more 
rapid . 

C. PECB1EM ENTEY 

Problem entry is accomplished by editing the main 
program directly. As an example, consider the following 
simple NIP with two design variables, and three constraint 
functions. 



2 

Minimize F(X)=X + 3XX + 2X 

1 12 



2 

- X - X + 1 

2 1 2 



subject to 



X + X - 3 < 0 
1 2 



1 



1 



+ 



2 < 0 



X 



X 



1 



2 
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2 < 0 



2 

X + X - X - 
1 1 2 

X >0.1 
i - 

With the MS COP loaded into memory and listed on the CRT, 
modifications are made on the program lines as follows to 
input this example : 

Line 100 

Just after the word "data", three integers are added, 
separated by a comma. The first number is 1JDV which is 
the number of design variables, the second is NIQC which 
is the number of inequality constraints, and the third is 
IPET which is print control number ( 0 ; only final 

results, 1 ; given data and final results, 2 ; given data 
and iterative subcptimal results) 

for example : 

100 data 2,3,2 

Lines 201-220 

Each line here corresponds to a separate design variable, 
beginning with X (1) and continuing in order to input 
X (NDV) . On each line, three values are separated by 
commas. After the word "data", these values are the 
initial values of the design variable, the lower bound on 
the variable and the upper bound on the variable. If no 
bound is to be specified, the entry is filled by "no". 

For the sample problem, the input is : 

201 data 3. , 0 . 1 , no 

202 data 3 . , 0 . 1, no 
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Lines 400 - 450 



These lines are available for defining the objective 
function. The objective function must be defined in 
terms of subscripted design variables X ( 1 ) , X (2) , etc. 

For the sample problem, the input is : 

400 f n_f = x (1) **2 + x (1) *x (2) + 2. *x (2) **2-x ( 1) -x (2) +1 . 

Lines 500-650 

These lines are available for defining the inequality 
constraint functions, which must be expressed using the 
format : 

601 if i = k then fn_g = G (x) - b. 

i i 

For the sample problem, the input is : 

00601 if i = 1 then fn_g = x(1)+x(2)-3. 

00602 if i = 2 then fn g = 1 . /x ( 1)+ 1 ./x (2) -2 . 

00603 if i = 3 then fn^g = x ( 1) **2 + x (1 ) -x (2) -2. • 

If there are many constant values in the constraint func- 
tions, the user may input data for these functions on 
lines 501-600 in order to simplify their statements. 
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IV. EXAMPLE PROBLEMS 

A. DESIGN OF CANTILEVERED BEAM 

1 * Dnif or m Cantilevered Beam 

Assume a cantilevered beam as shown in Figure 4. 1 
must be designed. The objective is to find the minimum 




r B 

i — n 

H 

l_i 



ij = 20000 Psi 
E = 30 E 6 Psi 
y = 1.0 inch 
P = 10000 lbs 




Figure 4.1 Design of a Uniform Cantilevered Beam. 



volume of material which will support the load P. 

The design variables are the width B and height H in 
the team. The design task is as follows : Find B and H to 



minimize volume 



V 



B H 1 



(4. 1) 



we wish to design the team subject to limit on bendin 
stress, shear stress, deflection and geometric conditions 
The bending stress in the beam must not exceed 20,000 psi . 

Me 6 P 1 

<T = = < 20,000 (4. 2) 

b I 2 

B H 



The shear stress must not exceed 10,000 psi. 




3 P 
2 A 



3 P 

< 10,000 

2 B H - 



(4. 3) 



and the deflection under the load must not exceed 1 inch. 



3 

s. -Li_ 

3 E I 



3 

4 P 1 



3 

E B H 



< 1.0 



(4.4) 



Additionally, geometric limits are imposed on the team size 
0.5 < B < 5.0 (4.5) 



1 .0 < H < 20.0 (4. 6) 

H/b < 1 0. (4. 7) 

Now we can input this problem to MSCOP. 

Input NDV , NIQC, IPRT 

00100 data 2,4,2 

Initial starting points 

00210 data 3. 5, 0.5, 5.0 
00220 data 16.0,1.0,20.0 

Evaluation of objective 

00400 f n_f = tl*x ( 1) *x (2) 
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Evaluation of constraints 



00500 


tl 




200 


• 




005 C 1 


be 


— 


30. 


e + 6 




00502 


be 


— 


10000. 




00503 


if 


i 


= 1 


then 


f n_g = 


00503 


if 


i 


= 2 


then 


fn_g = 


00503 


if 


i 


= 3 


then 


f n g = 


00503 


if 


i 


= 4 


then 


f n~g = 



6. *bp*tl/(20000 . *b*h*' 1 2) -1 . 

3. *bp/( 10000. *2. *b*h) - 1 . 

4. *bp*tl**3/ (be*b*h* *3) -1 . 
h/fc-30. 



TABLE I 

The Solution of a Uniform Cantilevered Beam 



objective ; 6664.0 



design variable ; 

X<1) = 1.852 
X (2) = 17.99 

constraint ; 



g { 1 ) = 0.000902 
g (2) = -0. 9549 
g (3) = -0.0109 
g (4) = -0.0286 



As a result of this problem are in Table 4.1. 

2. Variah le Cantilevered Beam 

The cantilevered beam shown in Figure 4.2 is to be 
designed for minimum material volume. The design variables 
are the width b and height h at each of 5 segments. We 

wish to design the beam subject to limits on 

stress (calcula te d at left end of each segment), deflection 
under the load, and the geometric requirement that the 
height of any segment does not exceed 20 times the width. 
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p 




y 




Cross section 



P = 50.000 N 
£ = 200 GHa 
L = 500 cm 
5 ~ 1 4 ,i x >0 N cm* 
v = 5 cm 



Figure 4.2 Design of a Variable Cantilevered Beam. 

The deflection y at the right end of segment i is 
calculated by the following recursion formulas : 




(4.8) 



y 



P 1. 

l 



E I 

i 



1 

i i 

L + + ZL 1 

2 j= 1 i 




(4.9) 



y 



2 

P 2 

i 



2 E I 

i 



2 1 

i i 

L - XT 1 + 

j=1 i 3 



+ y. i + 

1-1 1 



y i-1 



(4. 10) 
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where the deflection y is defined as positive downward, y’ 
is the derivative of y with respect to the X, and 1 i is the 
length of of segment i. Young’s modulus E is the same for 
all segments, and the moment of inertia for segment i is 



I '= 
i 



3 

t h 
i i 

12 



( 4 . 11 ) 



The tending moment at the left end of segment i is calcu- 
lated as 



fl = ? 

i 



L + 1 - H 1 

i j=1 i J 



(4. 12) 



and the corresponding maximum bending stress is 



C T = 

i 



M h 

1 i 

2 I 



(4. 13) 



The design task is now defined as 

N 



Minimize 



Subject to : 



V = ZL b h 1 
i=1 i i i 



(4. 14) 
(4. 15) 



0"c 



- 1 < 0 



i = 1 , . . . , N 



(4. 16) 



- 1 < 0 



(4. 17) 



h -20 b < 0 i = 1,...,S 

i i - 



(4. 18) 
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b >1.0 
i - 



h >5.0 

i - 



= 1 , N 



( 4 . 19 ) 



where <? is the allowable bending stress and y is the allo- 
wable displacement. This is a design problem in 10 vari- 
ables. There are 6 nonlinear constraints defined by Eg (4. 16) 
and Eg (4. 17), and 5 linear constraints defined by Eg (4. 18), 
and 10 side constraints on the design variables defined by 
Eg (4 . 19) . 

Now we can input this problem to NSCOP. 

Input NDV, NIQC, IPRT 

001 00 data 10,1 1,2 

Initial starting points 

00210 data 5.,1.,no 
00220 data 5.,1.,no 
00230 data 5.,1.,no 
00240 data 5.,1.,no 
00250 data 5.,1.,no 
00260 data 40. ,5., no 
00270 data 40. ,5., no 
00280 data 40. ,5., no 
00290 data 40., 5., no 
00300 data 40. ,5., no 



Evaluation of objective 

00400 f n_f = 100. * ( x ( 1 ) * x (6) 

x (4) *x (9) + x (5) *x ( 1 0) ) 



+ x (2) * x (7) + x (3) *x (3) 



Evaluatiom of constraints. 



004 90 def fn g (x, i) 

00498 dim bm710) £ bi (10) , 
00500 pcfc = 50000. 



sigi ( 10) ,ypb (10) ,yb (10) 



00501 be = 200. e+5 

00502 tl = 200. 

00503 sigb = 14000. 

005C4 ytb = .5 

00505 si = 40. 

00506 fcr m = 1 to 5 

00507 bm (m) = pet* (tl+sl-m*sl) 

00508 next m 

00509 for m = 1 to 5 

005 1 0 km = m + 5 

005 1 1 bi(m) = x (m) * x (km) ** 3/1 2 . 

C0512 sigi (m) = bm (m) *x (km) / (2. *bi (m) ) 

00513 next m 

00514 yzo = 0 . 

00515 ypzo = 0. 

00516 for m = 1 to 5 
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00517 

005 18 

005 1 9 

00520 

00521 

00522 

00550 

00560 

00570 

00580 

00590 

00600 

006 1 0 

00620 

00630 

00632 

00634 

00636 



ypb(m) = ( pcb*sl* (tl+sl/2 . -m*sl) ) / (be* bi (m) ) -ypzo 

aumi = pcb *£l**2* ( t 1-ra* si + 2. *sl/3 . ) 
yb (m) = dumm/ (2. *be*fci (m) ) +y pzo*sl+yzo 

ypzo - ypb (ir) 
yzo = yb(ra) 
next m 

rem constraint function, 
if i = 1 then fn_g = sigi Ml /sigb-1 . 



if i = 2 then fn_g 

if i = 3 then fn_g 

if i = 4 then fn_g 

if i = 5 then fn_g 

if i = 6 then fn_g = 

if i = 7 then fn_g = 

if i = 8 then fn_g = 

if i = 9 then fn g = 

if i = 10 then fn_g = 

if i = 11 then fn_g = 




2 \ /sigb- 1 . 

3) /sigb- 1 . 

4) /sigb- 1 . 



= si gi (5) /sigb- 1 



YM5)/ X bb-1 



*x(1] 

7) -20 . *x 2 

8) -20 . *x (3 
(9) -20. *x (1 
(ID) -20.*x (5) 






TABLE II 

The Solution of a Variable Cantilevered Beam 



objective ; 62133.35 

design variables constraints 



X{1) 


•= 


2.994 


G ( 1) 


= 


-0. 00219 


X(2) 


= 


2.782 


G (2) 


= 


-0. 00415 


X{3) 


- 


2. 528 


G (3) 


= 


-0. 00508 


X(4) 


- 


2.203 


G (4) 




-0.00406 


X (5) 




1.761 


G (5) 




-0.0177 


X<6) 




59.88 


G(6) 


= 


-0. 440 1 


X(7) 


- 


55.62 


G(7) 


= 


-0. 0 101 


X(8) 


- 


50. 56 


G (8) 


- 


-0. 0231 


X (9) 


= 


44. 14 


G<9) 




0. 0000 


X(10) 


= 


3 5. 19 


G { 1 0 ) 




-0. 0248 








G ( 1 1 ) 


- 


-0- 0278 
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E 



SIMPLE TRUSS 



Y 




P = 50000 N 
E = 200 GPa 

+ 14000 N/cm 
l = 100 cm 



Figure 4.3 Design of a 5-Bar Truss. 

A simple truss with 5 members as shown in Figure 4.3 is 
•designed for the minimum volume. The design variables are 
the sectional areas of the members. The constraints are 
formed for the stresses of the members not to exceed the 
given allowable stress. The lower bound for each design 
variable is also considered. The stresses are obtained by 
the displacement method of the finite element analysis. The 
equation to be solved is given by 

K-u = P (4.20) 

where K is the stiffness matrix, u. is the displacement 
vector and P is the lead vector as follows : 
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u 



(4. 21) 



u ] 




0 


1 






V 




0 


1 


P = 


0 


u 




2 






V 




-5000 


2 







A, L As A s 
1 VI JL <*- 



0 0 



K 



_ E 



As 
VI 5. 



0 



0 



Aa Aj 

I Vafi- 



0 



/1l 

X 



o 



A? A» _ /u 
JL r VIS. X 



Aa _ ^ _A* 

2- *. jl + vs: ft. 



From Eg. (4.20) the displacements are solved by 

-1 

0 = K -P 

Having displacements at all nodes, we can calc 
stress for each element. 



E* A1 



01= E •£ = 



where 



Al = 



A1 = 



I 2 2 

J( 1 + u ) + v 

V 1 1 1 

i ( v v T 



2 2 
) + ( u - u ) - 1 

2 12 2 



A1 



= i 



2 2 
(l+u) + v - 1 

3 2 2 3 



( 4 . 12 ) 



(4.23) 
ulate the 



(4. 24) 



(4.25) 
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/ 2 2 
A1 , = A 1 - t * + ( 1 o' V J " 1 

4 « 3 2 2 2 4 



2 2 
A1 =/(l+u) + ( 1 + v ) -1 

5 J 3 1 2 1 5 



The design problem is given by 



minimize V = ^ A 1 

i= 1 i i 



(4.26) 



Subject tc 

f cr, 



G = 
i 



- 1 . 0 < 0 i = 1 , .. . , 5 (4. 27) 



A >0.1 

i - 



i = 1 ,. . . , 5 (4.28) 



The MSCOF input for this problem is given as follows : 
Input NDV, NIQC, IPET 
00100 data 5, 5, 2 



Initial starting point 



00200 data 3.,.1,no 
00202 data 3.,.1,no 
00204 data 3 . , . 1 , n o 
00206 data 3.,.1,no 
OO208 data 3. , . 1, no 



Evaluation of objective 

00400 fn f = 100 * ( x(1) + x (2) + x(3) + sqr(2.)*x(4) + 

scr (2.) *x"(5) ) 



Evaluation of constraints 



0 50 0 dim vv (5) 

050 1 te = 2. e+ 7 

0502 tl = 100. 

0503 sigb = 14000. 
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0504 

0505 

0506 

0507 
0506 
0509 

05 10 
051 1 
0612 

0513 

0514 

0515 

0516 

0517 
0516 
0515 
0520 
0 52 1 

0522 

0523 
0690 
0592 

0594 

0595 

0596 
0598 
0600 
0602 
0604 
0606 
0608 
0610 

06 5 0 



cs = 2. *sqr (2 . ) 
ct = te/tl 

m : wnmv cs, ‘ ct 

k i 2 

+ x (5)/cs) *ct 
❖ ct 

+ x (4)/cs) *ct 
*ct/cs 



k 2 1 = 
k22 = 
k24 = 
k33 = 
k 34 = 
k42 = 
k43 = 
k 44 = 
d k 1 = 
dk2 = 
dk3 = 



X I Z 

111 



k2 4 
k34 



iin^i 4)/cs ’ 



*ct 



.. .1/k1_ 

- (k12+k22*dk1) /k24 
-k3 3 *dk 2 /k3 4 

pp = -60000 

v v 
vv 

vv (3) - - 

vv (4) = dk2*vv ( 

1 = s ?r ( (tl + vv 
tl + vv 
tl+ vv 
) *tl 



(1) = pp/(k42*c 
(21 = dk 1*vv (1) 
v (3) = dk3*vv 



2*dk1+k43*dk3 + k44*dk2) 



dll = sir 
dl2 = sqr 
d 13 = sqr 
hi = sqr ( 




**2) -tl 






- vv ( 3) ) ** 2) -tl 



dl4 = 
d 15 = 



sqr 

sqr 



l( 



hl + vv 
hl + vv 



reir constraint 



**2 + (hl-vv 
** 2+ (hl + vv 



i ** 2) -hi 
**2) -hi 



if 

if 

if 

if 

if 



f n end 



1 

2 

3 

3 

5 



t hen 
t hen 
t hen 
then 
then 



f n g 
fn_g 
f n g 
fn g 
f n~g 



= te*dl 1/ 
= te*dl2/ 
= te*dl3/ 
= t e * d 1 4 / 
= te*dl5/ 



tl*sigb) - 1 . 
tl*sigb) - 1 . 
4 tl*sigb) - 1 . 
/hl*sigb) -1. 
(hl*sigb) - 1. 



TABLE III 

The Sclution of a 5-Bar Truss 



objective ; 108-52 




design variables 


constraint 


X (1) = 0-1 


G (1) = - 1. 9988 


X (2) = 0.1 


G ( 2) = -2-0030 


X (3) = 3.514 


G (3) = -0.0030 


X ( 4) = 4.948 


G (4) - -0. 120 3 


X (5) = 0.1 


G ( 5) = - 1. 8797 
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V. SUMMARY AND CONCLUSION 



Numerical optimization is a powerful technique for those 

* 

confronted with practical engineering design problems. It 
is also, a useful tool for obtaining reasonable solutions to 
the classical engineering design problems. Since many engi- 
neers are now using nicrocomputers for solving design prob- 
lems, the development of microcomputer software which can be 
easily used is needed. 

In this thesis, an algorithm for constrained optimiza- 
tion problems is programmed in standard BASIC language 
(WBASIC version 2.0) on an IBM 3033. The users can easily 
convert this to other microcomputers. 

MSCOE (Microcomputer Software for Constrained 
Optimization Problems) employs the method of feasible direc- 
tions and specific modifications of a one-dimensional search 
for constrained optimization. MSCOP has been validated by 
tests on three constrained optimization problems. Its 
performance is good and could be made better through refine- 
ment of the algorithm. 

Since microcomputers are available with reasonable 
memory size and computational speed, their capabilities will 
continue to improve as more engineering software becomes 
available. MSCOP is considered to be a first step toward 
more widespread use cf optimization techniques on microcom- 
puters. 
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APPENDIX A 



00 10 
0020 
00 21 
0030 
0040 
0050 
0060 
0070 
0030 
0090 
0100 
0110 
01 1 5 
0120 
0125 
0130 
0 135 
0140 

0150 

0160 
0200 
0210 
0360 
0370 
0375 
0380 
0390 
0400 
0410 
0420 
0430 
0440 
0450 
0480 
0490 
0500 
0510 
0520 
053 0 

0540 

0550 

0560 

06 5 0 
0700 

07 10 
0720 
0730 
0740 
0750 
0760 
0770 
0780 
0785 
0790 
0800 



MSCOP PBOGEAM LISTING 



( 6 ) 



option base 1 

dim x (21) ,x0 (211 , gey (5 11 ,ngcv (51) ,df (21) ,dg(51,21) 
dim t hta (2 1) , wrky (5 1. 5 1) 

dim a (51,2 V ,b (51,51;, , p (21 ) , y ( 2 1) , s f 2 1) , u ( 51 ) , c (5 1) 
dim lwrk (51) , n vrk (51) ,wrk1 (51) , wr*2 (ol) ,wrk3(5l) 
dim wrku(51),wrkl(51) ,lowb^21) ,uprb(21) , lo3 (6) , up? i 
rem input data 
gosub 10000 

rem input number of design variables and constraints 

read nav , n igc , iprt 

data 2,4,2 

for i = i to ndv 

rem input initial value of design variables 
read x (i) 
xO (i) = x (i) 
if nige = 0 then 160 
read Io$,up$ 

if loS = ’no' then lowb (i) = bnlo else lowb(i) = 

value (lo$) 

if upS = ’no’ then uprb(i) = bnup else uprb (i) = 
value (up$) 
next i 

d ata 3.5, 0.5, 10. 
data 16. ,1.0,20. 

rem evalute the objective-function 
obj = fn f (x) 
itri = 1“ 

rem objective function 
def fn f (xl 

fn f = 200. *x (1) *x (2) 
fnend - 

rem evaluate the constraints 
for i = 1 to niac 

qcv (i) = f n_g (x,i) 

next i 

rem constraint functions 
(X/i) 

be 



ir i 



2U0. 

30. e + 6 

10000 . 

= 1 then 



if 

if 



= 2 
= 3 



then 

then 



fn_g = 

fn g = 
fiTg = 



6. *bp*tl) / 

20 00 0.*x(1)*x (2) **2) -1 . 

3. *bp) / (20000. *x (1) *x (2) ) -1 

4. *bp*tl**3) / 

be*x ( 1) *x (2) **3) -1. 
x (2) /(10 . *x ( 1) ) - 1. 



if i = 4 then fn g = 
fnend 

rem initial counting number input 
ical = 1 

if ical > 3 then stop 

rem call the optimization code. 

gosub 2000 

rem print results. 

rem 

rem re-counting number input, 
ical = ical+1 
if ical = 3 then 850 

rem 10% reduce the design variables, 
for i = 1 to ndv 



0810 

0320 

0830 

0840 

0850 

0860 

0870 

0880 

0890 

0900 

2000 

200 1 

2002 

2003 

2004 

2003 

2010 

2020 

2030 

2040 

2050 

2060 

2070 

2080 

2090 

2100 

2110 

2120 

2130 

2140 

2150 

2160 

2170 

2180 

2190 

2200 

2210 

2220 

2230 

224 0 

2250 

2260 

2270 

2280 

229 0 

2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2450 

2460 

2470 

2480 

2490 

2500 

3000 



x (i) = 0.9* 
xO (i) = x (i 
next i 
goto 720 
rem 10T increa 
for i = 1 to n 
x (i) = 1.1* 
xO (i) = x (i 

next i 
goto 720 
rem calculate 
obi = fn f (x) 
for i = 7 to n 
gcv(i) = fn 
next i 
itrg = 1 
itrg = itrg+1 
rem calculate 
constraint 
gosub 3500 
rem calculate 
active or 
gosub 3800 
if nave = 0 th 
gosub 3900 
rem calculate 
gosub 4000 
rem making the 
rem normalized 
gosub 4100 
rem normalized 
gosub 4200 
if nvc > 0 the 
rem calaulate 
gosub 5000 
goto 2230 
rem normalize 
for i = 1 to n 



x (i) 



s € design variables, 
dv 
x (i) 



the obj. constraint fen. 
igc 

_9 (*,i) 



the number of active and viola te 

s. 

the gradient of objective and 
violated constraints. 

en 2190 

the push-off factors 

matrix c 
the df (i) 

the DG (i) 

n gosub 4400 else gosub 4600 

the usable-feasible direction s (i) 



the df (i) 
dv 



s(i) = -(df(i).J 
next l 

rem normalize the s (i) 
gosub 5700 

rem one-dimensicnal search 

if nvc = 0 then gosub 6000 else gosub 9000 
rem update x for alph 
gosub 7000 
gosub 7100 

rem calculate new point value, 
not j = fn_f (x) 
rem convergence test 
gosub 6780 

if walp <= accx and delf <= dabf then 2470 
itn = itri + 1 

if itri > mxit then print ’check the problem’ 
obj = nobj 
for i = 1 to ndv 
x0.(i) = x (i ) 
next l 

for i = 1 to nige 

qcv(i) = fn_g(x,i) 
next i 

if iprt = 2 then 2460 
gosub 9200 
goto 2010 

rem print final results 
print ****** final results ***** ’ 
gosub 9200 
return 

rem initialize the integer working array 
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3005 
3010 
3015 
3020 
3050 
3 05 5 
3060 
3065 
3070 
3100 
3105 
3110 
311 5 
3120 
3150 
3155 
3160 
3165 
3170 
3200 
3205 
3210 
3215 
3220 
3250 
3255 
3260 
3265 
3270 
3275 
3280 
3300 
3305 
3310 
3315 
3320 
3350 
3353 
3356 
3359 
3362 
3365 
3368 
337 1 
3374 
3377 
3330 
3383 
3400 
3405 
3410 
34 1 c 
3420 
3425 
3430 
3450 
3455 
3460 
3465 
3470 
3475 
3480 
350 0 

3502 

3504 

3510 

3520 

3530 



for i = 1 to niqm 
iwrk (i) = 0 
next i 
return 

rem initialize the integer working array 
for i = 1 to niqm 



jwrk (i) = 0 

rt i 



nexi 
return 

rein initialize, the one-dimension working array 
for i = 1 to nicm 
wrkl (i) = 0 . 
next i 
return 

rem initialize the one-dimension working array 
for i = 1 to niqm 
wrk2(i) = 0. 
next i 
return 

rem initialize the one-dimension working array 
for i = 1 to niqc 
wrk3 (i) = gcv (i) 
next i 
return 

rem initialize the two-dimension working array 
for i = 1 to niqm 
for j = 1 tc ndvm 
wrky (i / j ) = 0. 
next g 
next i 
return 

rem initialize the derivative of objective DF (i) 
for i = 1 to ndvm 
df (i) =0. 
next i 
return 

rem initialize the a (i , j) , p (i) , y (i) , c (i) 
for i = 1 to ndvm 
P(i) = 0. 

7 ( 1 ) = 0 . 

:or j = 1 tc niqm 
a (j,i) = 0. 

next j 
next i 

for j = 1 to niqm 
c( j) = 0. 
next j 
return 

rem initialize the derivative of constraints DG(i,j) 
for i = 1 to niqm 
for j = 1 to ndvm 
dg (i, j) = 0. 

next j 
next i 
return 

rem initialize the b(i,j) 
for i = 1 to niqm 
for j = 1 to'niqm 
b (i, j) =0. 
next j 
next i 
return 

rem Calculate the number of active and violate 
constraints, 
gosub 3000 
gosub 3100 
nac = 0 
nvc = 0 

for i = 1 to niqc 
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3540 
3550 
3560 
3570 
3530 
3590 
36 10 
3620 
3630 
3640 
3650 
3660 
3670 
3680 
3690 
3700 
3710 
3720 

3 73 0 
3740 
3750 
3790 
3800 
3805 
3810 
3815 
3820 
3825 
3830 
383 5 
3840 
3850 
3860 
3900 
390 1 
3910 
3915 
3920 
3925 
3930 
3935 
3940 
3945 
3950 
3955 
3960 
3966 
4000 
4010 
4020 
4030 
4040 
4090 
4100 

4 102 
4105 
4 110 
4115 
4120 
4125 
4127 
4130 
4135 
4 140 
4145 
4200 
4202 
4205 
4210 



if gcv jij 



>= vcc then 3580 
if gcv (i) < acc then 3590 

nac = nac +1 
goto 3590 
nvc = nvc +1 
next i 

nave = nac+nvc 
if nave = 0 then 3790 
ii = 1 
1 j '= 1 

for i = 1 tc nige 

if gcv(i) >= vcc then 3720 
if gcv (i) < acc then 3750 
iwrk(nvc+ii) = i 
wrkl (nvc+ii) = gcv(i) 
ii = ii+ 1 
goto 3750 
lwrk ( j j) = i 
wrkl (gg) = gcv (i) 

3D.= 33 + 1 
next i 
return 

rem calculate the gradient of f (x) 

gosub 3300 

for i = 1 to ndv 

dxi = fdm*abs (x (i) ) 

if dxi <= mfds then dxi = mfds 
x (i) = x (i) +dxi 
dobg* = fn f (x) 
df (l) = (dobi-ob j) /dxi 
x (i) = xd (i) 
next i 
return 

rem calculate the DG(i,j) 
gosub 3400 
for i = 1 to ndv 
dxi = f dm*x (i) 

if dxi < mfds then dxi = mfds 
x (i) = x fi) +dxi 
for j = 1 tc nave 
k = iwrk ( 3 ) 
dcon = fn_g(x,k) 
dg(j,i) = (aeon- wrk 1 (j) ) /dxi 
next g 

x (i) = xO (i) 
next i 
return 

rem calcilate the push-off factor 
for i = 1 to nave 

thta(i) = t htO* (1 . -wrkl (i) /acc) 
if thta (i) > thtm then thta (i) 

next i 
return 

rem normalize the DF(i) 
gosub 3200 
fsg = 0. 

for i = 1 to ndv 

fsg = fsg + df(i )**2 
next i 

fsg = sgr(fsg) 

if fsg = 0 . then fsg = zro 

for i = 1 to ndv 

wrk3 (i) = ( 1 ./fsg) *df (i) 
next i 
return 

rem normalize the DG (i) 
gosub 3250 
for i = 1 to nave 
gsg = 0 . 



= thtm 
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4215 

4220 

4225 

4230 

4232 

4235 

4240 

4245 

4250 

4255 

4400 

4405 

4410 

4420 

4430 

4440 

4450 

4460 

4470 

4480 

4490 

4500 

4510 

4520 

4530 

4540 

4550 

4560 

4570 

4580 

4536 

4590 

4600 

4605 

46 10 

4620 

46 3 0 

4640 

4650 

46 6 0 

4670 

4680 

4690 

4700 

4710 

4720 

4730 

4740 

4750 

4760 

4770 

5000 

5002 

5005 

5010 

5040 

5050 

5060 

5 07 0 

5080 

5090 

5100 

5110 

5120 

5130 

5140 

5150 

5160 

5170 



for j = 1 tc ndv, 



nex 



cjsg = gsq+dg (i, j) 



gsg = sgr (gsg) 
if gsg = 0. tm 



1 = 1 tc ndv 
ai/j) = wrky(i,j) 



fcir 



gsg = u. tnen gsg = zro 
for j = 1 to ndv 

wrky (i, j) = ( 1. /gsg) *dg (i, j) 
next j 
next i 
return 

rent exist the violate constraints 
gosub 3350 
for i = 1 to nave 
for 

a 

next 1 

a(i,nav + 1) = thta(i) 

next i 

for i = 1 to ndv 
P (i) = -wrk3 (i) 
next i' 

(ndv+1) = phid 
■■ i = 1 to nave 

j = 1 tc ndv+1 
xx = a (i, j) *p (j) 
yy _= yy+xx 
next i 

c(i) = (-i.)*yy 

next i 
ndt = nave 
return 

rem only exist active constraints 
gosub 3350 
for i = 1 to nave 
for i = 1 tc ndv 

a(i/j) = wrky (i , j) 
next i 

a (i, ndv + 1) = thta(i) 

next i 

for j = 1 to ndv 

a(navc+1,j) = wrk3(j) 
next j 

a (navc+ 1 , ndv+ 1 ) = 1. 

p (ndv + 1) = 1 . 

for i = 1 to navc+1 

cc = a (i, ndv + 1 ) *p (ndv + 1) 
c(i) = (-1-)*cc 
next l 

ndt = navc+1 
return 

rem calculate the usable-feasible direction 

gosub 3000 

gosub 3250 

gosub 3450 

for i = 1 to ndt 

for j = 1 to ndv+1 
wrky ( j/ i ) = a (i, j) 
next j 
next i 

for i = 1 to ndh 



:or 



to ndb 
1 to ndv+1 



r 1 = 1 
ff = 0. 
for k = 

tf = a (i . k) * wrky ( k, i ) 
f f = ff+tf 
next k 

b(i,j) = (~1.)*ff 

next j 
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(-1.) *f f 



5130 
5190 
5200 
5210 
5220 
5230 
5240 
5250 
5260 
5270 
5280 
5290 
5300 
53 10 
5320 
5330 
5340 

53 5 0 
5360 
5370 
5380 
5385 
5400 
5405 

54 10 
5420 
5430 
544 0 
5460 
5470 
5480 
5490 
5500 
5505 
5510 
5520 
5525 
5530 

554 0 

555 0 
5560 
5565 
5570 
5580 
5600 
5610 
5620 
5630 

56 4 0 
5650 
5660 
5670 
5680 
5690 
5700 

57 10 
5720 
5730 
5740 
5750 
5755 
5760 
5770 
5780 
5820 
6000 
6005 
6010 
6015 



next i 

iter = 0 

nmax = 5*ndb 

rein begin iteration 

iter = iter+1 

cbmx = 0. 

ichk = 0 

for i = 1 to ndb 
ci = c ( i) 
bii = b (i , i ) 
if bii = 0. then 5340 
if ci > 0. then 5340 
cb = ci/bii 

if cb <= cbmx then 5340 



ichk 
cbmx 
next i 
if cbmx 
ichk 
33 = 

2i 



= l 

= cb 



if 



lr 



> nmax then 5550 



< zro or iter 
= 0 then 5550 
= iw rk (ichk) 

= 0 then iwrk (ichk) = ichk else iwrk(Lchk) = 0 
b (ichk, ichk) = 0. then b (ichk, ichk) = zro 
bb = 1 . /b (ichk, ich k) 
if bb = 0. then bb = zro 
for i = 1 to ndb 

b (ichk,i) = bb*b(ichk,i) 
next i 

c (ichk) = cbmx 
for i = 1 to ndb 

if i = ichk then 5530 
bbi = b (i , ichk) 
b (i , ichk) = 0. 
for j = 1 to ndb 
if i = ichk then 5520 

5 (i, j) = b (i, j) -bbi*b (ichk, j) 
n ext i 

c (i) = c (i) -bbi*cbmx 
next i 
goto 5220 
ner = 0 

for i = 1 to ndb 
u (i) = 0. 
j = iwrk(i) 

if 1 > 0 then u(i) = c(j) 
next i 

for i = 1 to ndb 
ff = 0. 

for i = 1 tc ndb 

fr = ff + vrky (i, j) * u ( j) 
next j 

ill) : ?!ir ££ 

next i 
return 

reir normalized the s(i) 

SSq = 0. 

for i = 1 to ndv 

ssq = ssq + s (i) ** 2 
next l 

ssq = sqr (ssq) 

if fslp = 0. then fslp = zro 
for i = 1 to ndv 

s (i) = (l./ssq) *s (i) 
next i 
return 

rem one-dimensicnal search 



rem calculate for slope oi 
fslp =0. 
for i = 1 to ndv 



for 
f (x) 



initial feasible ocint 
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6020 
6025 
6035 
604 0 
6045 
6050 
6055 
6060 
6065 
6067 
6070 

6075 

6076 
6030 
6085 
6090 
6095 
6100 
6105 
6110 
6115 
6120 
6125 
6130 
6135 
6140 
6145 
6150 
6155 
6160 
6165 
6170 
617 ^ 
6180 
6185 
6190 
6 19 5 

6200 
6205 
62 10 
6215 
6220 
6225 
6230 
62 3 5 
6237 
6240 
6245 
6250 
6252 
6255 
6260 
6265 
6270 
6275 
6280 
6285 
6290 
6295 
6300 
6305 
6310 
6315 

6320 
632 1 

6325 

6326 



fslp = fslp + df (i) *s (i) 
next i 

rem idenfy the initial point, 
alov = 0. 
flow = obj 
for i = 1 to nice 
wrkl (i) = gev (i) 
next i 

rem find alst : the 1st mid-point, 
if fslp = 0. then fslp = zro 
alst = aboj*f low/abs (fslp) 
for i = 1 to ndv 

if s(i) = 0. then s (i) = zro 
walp = alpx*x (i) /abs JsJi) ) 
if walp > alst then 60y5 
alst = walp 

next i 

rem update x for alst. 
alph = alst 
gosub 7000 
gosub 7100 

rem calculate the fist and wrkl(i) 
f 1st = f n_f (x) 
for i = 1 to nice 

wrkl (i) = f n_g (x,i) 
next i 

rem check the feasibility, 
ncvl = 0 

for i = 1 to niqc 

if wrkl (i) < vcc then 6170 

ncvl = ncvl+1 

next i 

if ncvl = 0 then 6200 
alst = 0 . 5*a Is t 
goto 6105 

rem find a2nd ; the 2nd mid-point, 
rem 2-points guadratic fit interpolation 
for minimum f (alpa) . 
aO = flow 
al = fslp 

a2 = (f 1st-a1*a1st-a0) /(a1st**2) 
if a2 <= 0. then a2 = zro 
a2nd = -al/ (2.*a2) 

rem 2-points linear interpolation for g(alpa)=0. 
for i = 1 to nice 
aO = wrkl (i) 

if alst = 0. then alst = zro 



al = (wrkl (i) -aO) /a 1st 
if al <= 0. then a 1 = 



zro 

walp = -aO/al 

if walp <= 0. then walp = 1000. 
if walp >= a2nd then 6265 
a2na = walp 

next i 

rem update x for a2nd. 
alph = a2nd 
gosub 7000 
gosub 7100 

rem calculate f2nd and wrk2 (i) 
f 2nd = f n_f (x) 
for i = 1 to niac 

wrk2 (i) = f n_g (x, i) 
next i 

rem find final point aupr by using 
3-points guadratic fit. 
fl = flow 
alpl = alow 
f 2 = fist 
alp2 = alst 



49 



6330 
633 1 
6335 
6340 
6342 
6345 
6347 
6350 
6355 
6360 
6 36 5 
63 7 0 
6 37 5 

6376 

6377 
6380 
6 38 5 
6390 
6395 
6400 
6405 
6410 
6415 
6420 
6425 
6430 
6435 
6440 
6445 
6450 
6455 
6460 
6465 
6470 
6475 
6480 
6485 
6490 
6495 
6500 
6505 
6510 

65 15 
6520 
6525 
6530 
6535 
6540 
6545 
6550 
6555 
6560 
6565 
6567 
6569 
6571 
6573 
6575 
6577 
6579 
6600 
6603 

6605 

6610 

66 15 
6620 
6630 



then 6380 
alpS 

fcr aupr 



fn_g (x,i) 
new point. 



£3 = f 2nd 
alp3 = a2nd 

rem 3-points quadratic fit interpolation, 
gosub 5600 

if a2 = 0. then a2 = zro 
a3rd = -a1/f2.*a2) 
if a3rd <= 0. then a3rd = 1000. 
for i = 1 to nigc 
f 1 = wrkl (i* 
f 2 = wrkl (i 
f3 = wrk2 (i 
gosub 6600 
gosub 6630 
if alps > a3rd 
a3rd = 

next i 

rem update x 
alph = a3rd 
gosub 7000 
gosub 7100 

rem calculate the fupr and wrku(i) 
fupr = £n_f(x) 
for i = 1 to niqc 
wrku(i) = 
n ext i 

rem find 4th 
fl = fist 
f 2 = f 2nd 
f3 = f3rd 
alpl = alst 
alp2 = a2nd 
alp3 = a 3rd 

rem 3-points quadratic fit. 
gosub 5600 

if a2 = 0. then a2 = zro 
aupr = -a1/(2.*a2) 
for i = 1 to niqc 
fl = wrkl (i' 
f 2 = wrk2 i' 
f 3 = wrk3 \i) 
alpl = alst 
alp2 = a2nd 
alp3 = a3rd 
gosub 6600 
gosub 6630 

if alps > aupr then 6540 
aupr = alps 
next l 

rem update x for aupr 
alph = aupr 
gosub 7000 
gosub 7100 

rem evaluate fupr and wrku(i) 
fupr = f n_f (x) 
for i = 1 to niqc 

wrku (i) = fn_g(x,i) 
next i 

rem find optimum alpa for not violating constraints. 

gosub 14300 

return 

rem 3-points quadratic fit. 
if alpl = alp2 cr alp2 = alp3 
then return 

a2 = ( (f 3-f 1 ) / (alp3-alp1 ) - 

(f2-f 1} / (alp 2-alp 1) / (alp3-alp2) 
a 1 = ]f2-f 1)/(alp2-alp1) -a2*(alp1+alp2) 
aO = f 1-a 1 *alp 1-a2*alp 1**2 
return 

rem zero of polynomial fcr g(alpa) 



or alpl = alp3 
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6635 

6640 

6642 

6645 

6650 

6655 

6660 

6665 

6670 

6675 

6680 

6685 

6690 

6695 

6700 

6705 

6710 

6712 

67 15 

6720 

6780 

6790 

6795 

6800 

6315 

6816 

6820 

6830 

6850 

6855 

6860 

6870 

6880 

6890 

6910 

6990 

7000 

7010 

7020 

7030 

7040 

7100 

7110 

7120 

7130 

7140 

7 150 

8000 

8010 

8020 

8030 

8040 

8050 

8060 

8070 

8080 

8090 

8100 

8110 

8120 

8130 

8140 

8150 

8160 

8170 

8180 

8190 

8200 

8210 



dd = a1**2-4.*a2*a0 
if dd < 0. then 6715 
if a2 <= 0. then a2 = zro 
if a2 = 0. then a2 = zro 

alpb = f-al +sgr (dd) ) / (2. *a2] 
alpc = (-a 1 -sqr jdd) ) / (2. *a 2 \ 

if alpb <= 0 ana alpc <= 0. 
if alpb >= 0. and alpc >= 0 
if alpb >= 0. and alpc < 0. 
alps = alpc 
goto 6720 

alps = alpb 
goto 6720 

if alpb >= alpc then 6710 
alps = alpb 
goto 6720 

alps = alpc 
goto 6720 
alps = 1000. 
return 

rein update aboj and alpx 

delf = abs (obi-nob j) 

aiff = abs (delf/objp 

abcj = fabo j+dif f ) /2. 

walp = 0. 

we lx = 0. 

for i = 1 to ndv 

delx = abs (xO (i) -x (i) ) 
difx = abs (delx/x0 (i) ) 
if delx >= welx then welx = 
if difx <= walp then 6880 
walp = difx 

next i 



alpx = (alpx+walp) /2. 
dabf = accf*abs (ob j) 
return 

rem update the x(i) 
for i = 1 to ndv 

x (i) = xO (i ) +alph*s (i) 
next i 
return 

rem check the side-constraints 
for i = 1 to ndv 

if x(i) < lcwb (i) then x (i) 
if x(i) > uprb (i) then x (i) 
next i 
return 

rem estimate the alpa 
fstr = flow 
alpa = alow 
nvcl = 0 

for i = 1 to nigc 

if wrkl (i) < vcc then 8070 

nvcl = nvcl+1 
next i 

if nvcl > 0 then 8120 
if fist > fstr then 8120 

alpa = alst 

fstr = fist 

nvcl = 0 

for i = 1 to nigc 

if wrk2(i) < vcc then 8160 
nvcl = nvcl+1 
next i 

if nvcl > 0 then 8210 
if f 2nd > fstr then 8210 

alpa = a2nd 

fstr = f2nd 

nvcl = 0 
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then 6715 
then 6695 
then 6685 



delx 



= lowb (i) 
= uprb (i) 



8220 
823 0 
8240 
8250 
82 6 0 
8270 
8280 
8290 
8300 
83 10 
8320 
8330 
8340 
8350 
8360 
83 7 0 
83 8 0 
3390 
8400 
9000 

9002 
9004 
9006 
9008 
9010 
9012 
9014 
9016 
9018 
90 20 
9022 
9024 

9026 

9027 

9028 
9030 
9032 
9034 
9036 
9038 
9040 
9042 
9044 
9046 
9048 
9050 
9052 
9054 
9056 
9058 
9060 
9062 
9064 
9066 
9068 
9070 
9072 
9074 
9076 
9078 
9030 
9200 
9205 
9210 
9215 
9220 
9225 
9230 



for i = 1 to nice 

if wrk3 (i) < vcc then 3250 

nvcl = nvcl +1 
next i 

if nvcl > 0 then 8300 
if f3rd > fstr then 8300 

alpa = a3rd 

fstr = f 3rd 

nvcl = 0 

for i = 1 to nige 

if wrku (i) < vcc then 8340 

nvcl = nvcl+1 
next i 

if nvcl > 0 then 8390 
if fupr > fstr then 8390 

alpa = aupr 

fstr = fupr 

aiph = alpa 
return 

rem one-di mens icnal search for initial 
infeasible point, 
ii = 1 

qcvm = wrkl ( 1) 
xor i = 1 to nave 
if wrkl (i) <= gevm then 9014 

ii = i 

qcvm = wrkl (i) 
next i 

rem calculate the slope of badly violation. 

qslp = 0. 

for i = 1 to ndv 

qslp = gslp+dg (ii, i) *s (i) 
next i 

rem calculate the alph. 

if qslp = 0. then gslp = zro 

alph = -gcvm/qslp 

rem update X for alph. 

gosub 7000 

gosub 7100 

rem evalute the objective and constraint, 
obj = fn f(x) 
for i = 7 to nige 

qcv (i) = fn_g (x,i) 
next i 

rem calculate the NVC. 
gosub 3500 

if nvc = 0 then return 
rem update initial value, 
for i = 1 to ndv 
x0.(i) = x (i) 

next l 

rem calculate df (i) ,dg (i , j) and push-off factor, 
gosub 3800 
gosub 3900 
gosub 4000 

rem normalize the df (i) , dg (i, j) 
gosub 4100 
gosub 4200 

rem find the search direction, 
gosub 4400 
gosub 5000 
goto 9000 

rem print the results 
print ' ’ 

print ************ data ************ 
print ' ’ 

print ’The number of design variables = ’.ndv 

print ’The number of inequality constraints = ’,niqc 
print ’ ' 
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9235 

9240 

924 5 

9250 

9255 

9260 

9265 

9270 

9275 

9280 

9285 

9290 

9295 

9 30 0 

9305 

9310 

9315 

9500 

9510 

9520 

9530 

9540 

9550 

9560 

9570 

9580 

9590 

9600 

96 10 

9620 

9630 

9640 

9650 

9660 

9670 

9680 

9690 

9700 

9800 



print 'The objective value = ',obj 
print " 

rint ****** design variables *****' 
or i = 1 to ndv 

print 'x(';i;') = ',x(i) 
next i 
print ' ' 

print 'the number of active constraints = ' ; nac 
print ' ' 

print ’the number of vionate constraints = ' ;nvc 
print ' ' 

print ***** constraint value ***** 
print " 

for i = 1 to niac 

print 'g(';i;') = ' ;gcv (i) 
next i 
return 

rem default number 

mxit = 50 ! maximum iteration number 

fdm = .01 ! finite difference step 

mfds = .001 ! maximum absolute finite difference step 
vcc = .004 ! violated constraint criteria (thickness) 

acc = -. 1 ! active constraints criteria (thickness) 

thtO = 1. ! push-off factor multiplier (theta zero) 

thtm = 50. ! maximum value of push-off factor 

phid = 100000. ! weighting-factor used in direction 

when infeasible 

accf = .001 ! absolute convergence criteria 

accx = 0.001 ! absolute convergence criteria.* 

zro = .0001 ! defined zero 

espl = 
bnlo = 
bnup = 
dalp = 

abcj = 
alpx = 
ndvm = 
nigm = 

return 
end 



.005 ! used to prevent division by zero 

-l.e+70 ! the value of low boundary 
l.e+70 ! the value of upper boundary 

.01 ! steD size of alpa m one-dimensional 

search 

0.1 ! step size for reduce objective 

.1 ! reduce the design variable factor 

21 ! the number of maximum design variable 

51 ! the number of maximum inequality 

constraints 
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